step2Model2ComputingEff_2016 <- function(sims){
  ATEm20<-rep(NA, sims)
  AMEm20<-rep(NA, sims)
  ADEm20<-rep(NA, sims)
  ATEm15<-rep(NA, sims)
  AMEm15<-rep(NA, sims)
  ADEm15<-rep(NA, sims)
  ATEm10<-rep(NA, sims)
  AMEm10<-rep(NA, sims)
  ADEm10<-rep(NA, sims)
  ATEm5<-rep(NA, sims)
  AMEm5<-rep(NA, sims)
  ADEm5<-rep(NA, sims)
  ATE0<-rep(NA, sims)
  AME0<-rep(NA, sims)
  ADE0<-rep(NA, sims)
  ATE5<-rep(NA, sims)
  AME5<-rep(NA, sims)
  ADE5<-rep(NA, sims)
  ATE10<-rep(NA, sims)
  AME10<-rep(NA, sims)
  ADE10<-rep(NA, sims)
  ATE15<-rep(NA, sims)
  AME15<-rep(NA, sims)
  ADE15<-rep(NA, sims)
  ATE20<-rep(NA, sims)
  AME20<-rep(NA, sims)
  ADE20<-rep(NA, sims)
  ATE25<-rep(NA, sims)
  AME25<-rep(NA, sims)
  ADE25<-rep(NA, sims)
  
  for (i in 1:sims){
    ## TOTAL EFFECT ##
    newxm25 <- cbind(t(holder[1:5,,(i-1)*11+1]),
                     rep(-2.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm20 <- cbind(t(holder[1:5,,(i-1)*11+2]),
                     rep(-2.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm15 <- cbind(t(holder[1:5,,(i-1)*11+3]),
                     rep(-1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10 <- cbind(t(holder[1:5,,(i-1)*11+4]),
                     rep(-1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5<- cbind(t(holder[1:5,,(i-1)*11+5]),
                   rep(-.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0 <- cbind(t(holder[1:5,,(i-1)*11+6]),
                   rep(0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5 <- cbind(t(holder[1:5,,(i-1)*11+7]),
                   rep(.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10 <- cbind(t(holder[1:5,,(i-1)*11+8]),
                    rep(1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15 <- cbind(t(holder[1:5,,(i-1)*11+9]),
                    rep(1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx20 <- cbind(t(holder[1:5,,(i-1)*11+10]),
                    rep(2.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx25 <- cbind(t(holder[1:5,,(i-1)*11+11]),
                    rep(2.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    ## MEDIATION EFFECT ##
    newxm25M <- cbind(t(holder[1:5,,(i-1)*11+1]),
                      rep(-1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm20M <- cbind(t(holder[1:5,,(i-1)*11+2]),
                      rep(-1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm15M <- cbind(t(holder[1:5,,(i-1)*11+3]), 
                      rep(-.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10M <- cbind(t(holder[1:5,,(i-1)*11+4]), 
                      rep(0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5M <- cbind(t(holder[1:5,,(i-1)*11+5]), 
                     rep(.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0M <- cbind(t(holder[1:5,,(i-1)*11+6]), 
                    rep(1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5M <- cbind(t(holder[1:5,,(i-1)*11+7]), 
                    rep(1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10M <- cbind(t(holder[1:5,,(i-1)*11+8]), 
                     rep(2.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15M <- cbind(t(holder[1:5,,(i-1)*11+9]), 
                     rep(2.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    ## DIRECT EFFECT ##
    newxm15D <- cbind(t(holder[1:5,,(i-1)*11+1]),
                      rep(-1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10D <- cbind(t(holder[1:5,,(i-1)*11+2]),
                      rep(-1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5D <- cbind(t(holder[1:5,,(i-1)*11+3]), 
                     rep(-.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0D <- cbind(t(holder[1:5,,(i-1)*11+4]), 
                    rep(0, nrow(df)),df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5D <- cbind(t(holder[1:5,,(i-1)*11+5]), 
                    rep(.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10D <- cbind(t(holder[1:5,,(i-1)*11+6]), 
                     rep(1.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15D <- cbind(t(holder[1:5,,(i-1)*11+7]), 
                     rep(1.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx20D <- cbind(t(holder[1:5,,(i-1)*11+8]), 
                     rep(2.0, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx25D <- cbind(t(holder[1:5,,(i-1)*11+9]), 
                     rep(2.5, nrow(df)), df$avgUnemp16, df$chgInc1615, df$inc16,  
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    
    
    coeffics <- matrix(as.numeric(maincoefs[i,]), ncol = 1)
    evalRandomState <- matrix(NA, nrow=length(groupcoefsID), ncol=2) # Create an empty matrix that will store random effects
    evalRandomState[,1] <- t(groupcoefs[i, ]) # Each row of groupcoefs has intercepts
    evalRandomState[,2]<-groupcoefsID[1:length(groupcoefsID)] # Col 2 is for stateID
    colnames(evalRandomState) <- c("Intercept", "stateID")
    newf<-merge(grouponly, evalRandomState, by = "stateID", all.x = TRUE)
    newf <- as.data.frame(cbind(newf$Intercept, newf$IDnum))
    names(newf) <- c("Intercept", "IDnum")
    newf <- newf[order(newf$IDnum),]
    
    etam25 <- newxm25%*%coeffics + newf[,1] # newf[,1]: random effects
    etam20 <- newxm20%*%coeffics + newf[,1]
    etam15 <- newxm15%*%coeffics + newf[,1]
    etam10 <- newxm10%*%coeffics + newf[,1]
    etam5 <- newxm5%*%coeffics + newf[,1]
    eta0 <- newx0%*%coeffics + newf[,1]
    eta5 <- newx5%*%coeffics + newf[,1]
    eta10 <- newx10%*%coeffics + newf[,1]
    eta15 <- newx15%*%coeffics + newf[,1]
    eta20 <- newx20%*%coeffics + newf[,1]
    eta25 <- newx25%*%coeffics + newf[,1]
    
    etaAMEm25 <- newxm25M%*%coeffics + newf[,1]
    etaAMEm20 <- newxm20M%*%coeffics + newf[,1]
    etaAMEm15 <- newxm15M%*%coeffics + newf[,1]
    etaAMEm10 <- newxm10M%*%coeffics + newf[,1]
    etaAMEm5 <- newxm5M%*%coeffics + newf[,1]
    etaAME0 <- newx0M%*%coeffics + newf[,1]
    etaAME5 <- newx5M%*%coeffics + newf[,1]
    etaAME10 <- newx10M%*%coeffics + newf[,1]
    etaAME15 <- newx15M%*%coeffics + newf[,1]
    
    etaADEm15 <- newxm15D%*%coeffics + newf[,1]
    etaADEm10 <- newxm10D%*%coeffics + newf[,1]
    etaADEm5 <- newxm5D%*%coeffics + newf[,1]
    etaADE0 <- newx0D%*%coeffics + newf[,1]
    etaADE5 <- newx5D%*%coeffics + newf[,1]
    etaADE10 <- newx10D%*%coeffics + newf[,1]
    etaADE15 <- newx15D%*%coeffics + newf[,1]
    etaADE20 <- newx20D%*%coeffics + newf[,1]
    etaADE25 <- newx25D%*%coeffics + newf[,1]
    
    
    ATEm15[i] <- mean(exp(etam15)/(1 + exp(etam15)), na.rm = T) - mean(exp(etam25)/(1 + exp(etam25)), na.rm = T)
    AMEm15[i] <- mean(exp(etam15)/(1 + exp(etam15)), na.rm = T) - mean(exp(etaAMEm25)/(1 + exp(etaAMEm25)), na.rm = T)
    ADEm15[i] <- mean(exp(etaADEm15)/(1 + exp(etaADEm15)), na.rm = T) - mean(exp(etam25)/(1 + exp(etam25)), na.rm = T)
    
    ATEm10[i] <- mean(exp(etam10)/(1 + exp(etam10)), na.rm = T) - mean(exp(etam20)/(1 + exp(etam20)), na.rm = T)
    AMEm10[i] <- mean(exp(etam10)/(1 + exp(etam10)), na.rm = T) - mean(exp(etaAMEm20)/(1 + exp(etaAMEm20)), na.rm = T)
    ADEm10[i] <- mean(exp(etaADEm10)/(1 + exp(etaADEm10)), na.rm = T) - mean(exp(etam20)/(1 + exp(etam20)), na.rm = T)
    
    ATEm5[i] <- mean(exp(etam5)/(1 + exp(etam5)), na.rm = T) - mean(exp(etam15)/(1 + exp(etam15)), na.rm = T)
    AMEm5[i] <- mean(exp(etam5)/(1 + exp(etam5)), na.rm = T) - mean(exp(etaAMEm15)/(1 + exp(etaAMEm15)), na.rm = T)
    ADEm5[i] <- mean(exp(etaADEm5)/(1 + exp(etaADEm5)), na.rm = T) - mean(exp(etam15)/(1 + exp(etam15)), na.rm = T)
    
    ATE0[i] <- mean(exp(eta0)/(1 + exp(eta0)), na.rm = T) - mean(exp(etam10)/(1 + exp(etam10)), na.rm = T)
    AME0[i] <- mean(exp(eta0)/(1 + exp(eta0)), na.rm = T) - mean(exp(etaAMEm10)/(1 + exp(etaAMEm10)), na.rm = T)
    ADE0[i] <- mean(exp(etaADE0)/(1 + exp(etaADE0)), na.rm = T) - mean(exp(etam10)/(1 + exp(etam10)), na.rm = T)
    
    ATE5[i] <- mean(exp(eta5)/(1 + exp(eta5)), na.rm = T) - mean(exp(etam5)/(1 + exp(etam5)), na.rm = T)
    AME5[i] <- mean(exp(eta5)/(1 + exp(eta5)), na.rm = T) - mean(exp(etaAMEm5)/(1 + exp(etaAMEm5)), na.rm = T)
    ADE5[i] <- mean(exp(etaADE5)/(1 + exp(etaADE5)), na.rm = T) - mean(exp(etam5)/(1 + exp(etam5)), na.rm = T)
    
    ATE10[i] <- mean(exp(eta10)/(1 + exp(eta10)), na.rm = T) - mean(exp(eta0)/(1 + exp(eta0)), na.rm = T)
    AME10[i] <- mean(exp(eta10)/(1 + exp(eta10)), na.rm = T) - mean(exp(etaAME0)/(1 + exp(etaAME0)), na.rm = T)
    ADE10[i] <- mean(exp(etaADE10)/(1 + exp(etaADE10)), na.rm = T) - mean(exp(eta0)/(1 + exp(eta0)), na.rm = T)
    
    ATE15[i] <- mean(exp(eta15)/(1 + exp(eta15)), na.rm = T) - mean(exp(eta5)/(1 + exp(eta5)), na.rm = T)
    AME15[i] <- mean(exp(eta15)/(1 + exp(eta15)), na.rm = T) - mean(exp(etaAME5)/(1 + exp(etaAME5)), na.rm = T)
    ADE15[i] <- mean(exp(etaADE15)/(1 + exp(etaADE15)), na.rm = T) - mean(exp(eta5)/(1 + exp(eta5)), na.rm = T)
    
    ATE20[i] <- mean(exp(eta20)/(1 + exp(eta20)), na.rm = T) - mean(exp(eta10)/(1 + exp(eta10)), na.rm = T)
    AME20[i] <- mean(exp(eta20)/(1 + exp(eta20)), na.rm = T) - mean(exp(etaAME10)/(1 + exp(etaAME10)), na.rm = T)
    ADE20[i] <- mean(exp(etaADE20)/(1 + exp(etaADE20)), na.rm = T) - mean(exp(eta10)/(1 + exp(eta10)), na.rm = T)
    
    ATE25[i] <- mean(exp(eta25)/(1 + exp(eta25)), na.rm = T) - mean(exp(eta15)/(1 + exp(eta15)), na.rm = T)
    AME25[i] <- mean(exp(eta25)/(1 + exp(eta25)), na.rm = T) - mean(exp(etaAME15)/(1 + exp(etaAME15)), na.rm = T)
    ADE25[i] <- mean(exp(etaADE25)/(1 + exp(etaADE25)), na.rm = T) - mean(exp(eta15)/(1 + exp(eta15)), na.rm = T)
    
    if (i %% 20 == 0) print(i)
    if (i %% 20 == 0) cat("Will be done when this number becomes", sims, ". Keep cool :) \n") 
  }
  return(as.data.frame(cbind(ATEm15, AMEm15, ADEm15, 
                             ATEm10, AMEm10, ADEm10, 
                             ATEm5, AMEm5, ADEm5, 
                             ATE0, AME0, ADE0, 
                             ATE5, AME5, ADE5, 
                             ATE10, AME10, ADE10,
                             ATE15, AME15, ADE15,
                             ATE20, AME20, ADE20,
                             ATE25, AME25, ADE25)))
}